import numpy as np
import matplotlib.pyplot as plt
#根据题目条件，确认常量为
m = 239.46
g = 9.8
V = 0.2058
r = 1035.71
k = 0.6
#使用四阶龙格-库塔方法求解微分初值问题
h = 0.002
def f(xn,yn):
    return g*(m-r*V-k*yn)/(m*yn)
    #return yn-2*xn/yn
def K1(xn,yn):
   return f(xn,yn)
def K2(xn,yn):
   return f(xn+h/2,yn+h/2*K1(xn,yn))
def K3(xn,yn):
   return f(xn+h/2,yn+h/2*K2(xn,yn))
def K4(xn,yn):
   return f(xn+h,yn+h*K3(xn,yn))
x = np.arange(0,100,h)
y = np.zeros(len(x))
y[0] = 0.001
for i in range(1,len(x)):
    y[i] = y[i-1] + h/6*(K1(x[i-1],y[i-1])+2*K2(x[i-1],y[i-1])+2*K3(x[i-1],y[i-1])+K4(x[i-1],y[i-1]))
    print(i,y[i],x[i])